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ABSTRACT 

The loss of linear momentum by gravitational radiation and the resulting gravita- 
tional recoil of black-hole binary systems may play an important role in the growth 
of massive black holes in early galaxies. We calculate the gravitational recoil of non- 
spinning black-hole binaries at the second post-Newtonian order (2PN) beyond the 
dominant effect, obtaining, for the first time, the 1.5PN correction term due to tails 
of waves and the next 2PN term. We find that the maximum value of the net recoil 
experienced by the binary due to the inspiral phase up to the innermost stable cir- 
cular orbit (ISCO) is of the order of 22kms~^. We then estimate the kick velocity 
accumulated during the plunge from the ISCO up to the horizon by integrating the 
momentum flux using the 2PN formula along a plunge geodesic of the Schwarzschild 
metric. We find that the contribution of the plunge dominates over that of the inspi- 
ral. For a mass ratio m2/mi = 1/8, we estimate a total recoil velocity (due to both 
adiabatic and plunge phases) of 100 it 20kms~^. For a ratio 0.38, the recoil is max- 
imum and we estimate it to be 250 it 50kms~^. In the limit of small mass ratio, we 
estimate V/c ^ 0.043 (1 it 20%) (m2/mi)^. Our estimates are consistent with, but span 
a substantially narrower range than, those of Favata et al. (2004) . 

^On leave from University of Sana, Yemen 

^Visiting Researcher at C/ReCO - Gravitation et Cosmologie, Institut d'Astrophysique de Paris, C.N.R.S., 98 bis 
boulevard Arago, 75014 Paris, France 
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1. Introduction and summary 

The gravitational recoil of a system in response to the anisotropic emission of gravitational 
waves is a phenomenon with potentially important astrophysical consequences (Merritt et al. 2004). 
Specifically, in models for massive black hole formation involving successive mergers from smaller 
black hole seeds, a recoil with a velocity sufficient to eject the system from the host galaxy or mini- 
halo would effectively terminate the process. Recoils could eject coalescing black holes from dwarf 
galaxies or globular clusters. Even in galaxies whose potential wells arc deep enough to confine 
the recoiling system, displacement of the system from the center could have important dynamical 
consequences for the galactic core. Consequently, it is important to have a robust estimate for the 
recoil velocity from inspiraling black hole binaries. 

Recently, Favata et al. (2004) estimated the kick velocity for inspirals of both non-spinning and 
spinning black holes. For example, for non-spinning holes, with a mass ratio of 1:8, they estimated 
kick velocities between 20 and 200kms~^. The result was obtained by (i) making an estimate of 
the kick velocity accumulated during the adiabatic inspiral of the system up to its innermost stable 
circular orbit (ISCO), calculated using black-hole perturbation theory (valid in the small mass ratio 
limit), extended to finite mass ratios using scaling results from the quadrupole approximation, and 
(ii) combining that with a crude estimate of the kick velocity accumulated during the plunge phase 
(from the ISCO up to the horizon) . The plunge contribution generally dominates the recoil, and is 
the most uncertain. 

It is the purpose of this paper to compute more precisely the gravitational recoil velocity 
during the inspiral phase up to the ISCO, and to attempt to narrow that uncertainty in the plunge 
contribution for non-spinning inspir ailing black holes. 

Earlier approaches for computing the recoil of general matter systems include a near-zone 

computation of the recoil in linearized gravity (Peres 1962), flux computations of the recoil as an 
interaction between the quadrupole and octupole moments (Bonnor & Rotenberg 1961; Papapetrou 
1962), a general multipole expansion for the linear momentum flux (Thorne 1980), and a radiation- 
reaction computation of the leading-order post-Newtonian recoil (Blanchet 1997). 

Using the post-Minkowskian and matching approach (Blanchet &: Damour 1986; Blanchet 
1995, 1998) for calculating equations of motion and gravitational radiation from compact binary 
systems in a post-Newtonian (PN) sequence, Blanchet et al. (2002, 2004) have derived the grav- 
itational energy loss and phase to 0[{vlcf\ beyond the lowest-order quadrupole approximation, 
corresponding to 3.5PN order, and the gravitational wave amplitude to 2.5PN order (Arun et al. 
2004). Using results from this program, we derive the linear momentum flux from compact binary 
inspiral to C'[(t;/c)'^], or 2PN order, beyond the lowest-order result. 

The leading, "Newtonian" contribution^ for binaries was first derived by Fitchett (1983), and 
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Fig. 1. — Recoil velocity as a function of rj. 



was extended to IPN order by Wiseman (1992). We extend these results by including both the 
1.5PN order contributions caused by gravitational-wave tail effects, and the next 2PN order terms. 
We find that the linear momentum loss for binary systems in circular orbits is given by ^ 
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where m = mi +m2, Sm = mi —m2, f] = mim2/m^ (we have < t] < 1/4, with rj = 1/4 for equal 
masses), and where x = (mw)^/^ is the FN parameter of the order of 0[{v/c)'^], where oj = dcji/dt is 
the orbital angular velocity. The quantity A* is a unit tangential vector directed in the same sense as 
the orbital velocity v = vi — V2. The term at order x^/^ = 0[(v/c)^] comes from gravitational- wave 
tails. Notice that, as expected for non-spinning systems, the flux vanishes for equal-mass systems 
{5m = or ?7 = 1/4). 



it really corresponds to a 3.5PN radiation-reaction effect in the local equations of motion. 

^In most of this paper we use units in which G = 1 = c. We generally do not indicate the neglected PN remainder 
terms (higher than 2PN). 
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To calculate the net recoil velocity, we integrate this flux along a sequence of adiabatic quasi- 
circular inspiral orbits up to the ISCO. Wc then connect that orbit to an unstable inspiral orbit of a 
test body with mass fi = rjm in the geometry of a Schwarzschild black hole of mass m, with initial 
conditions that include the effects of gravitational radiation damping. Using an integration variable 
that is regular all the way to the event horizon of the black hole, we integrate the momentum flux 
vector over the plunge orbit. Combining the adiabatic and plunge contributions, calculating the 
magnitude, and dividing by m gives the net recoil velocity. Figure 1 shows the results. Plotted as 
a function of the reduced mass parameter rj are curves showing the results correct to Newtonian 
order, to IPN order, to 1.5PN order and to 2PN order. Also shown is the contribution of the 
adiabatic part corresponding to the inspiral up to the ISCO (calculated to 2PN order). The "error 
bars" shown are an attempt to assess the accuracy of the result by including 2.5PN and 3PN terms 
with numerical coefficients that are allowed to range over values between —10 and 10. 

We note that the IPN result is smaller than the Newtonian result because of the rather large 
negative coefficient seen in Eq. (1). On the other hand, the tail term at 1.5PN order plays a 
crucial role in increasing the magnitude of the effect (both for the adiabatic and plunge phases), 
and wc observe that the small 2PN coefficient in Eq. (1) leads to the very small difference between 
the 1.5PN and 2PN curves in Fig. 1. In our opinion this constitutes a good indication of the 
"convergence" of the result. The momentum flux vanishes for the equal-mass case, r/ = 1/4, and 
reaches a maximum around 7] = 0.2 (a mass ratio of 0.38), which corresponds to the maximum of 
the overall factor rf5m/m = rf'^/l — 4r/, reflecting the relatively weak dependence on r] in the PN 
corrections. Wc propose in Eq. (36) below a phenomenological analytic formula which embodies 
this weak r] dependence, and flts our 2PN curve remarkably well. 

In contrast to the range 20 - 200 kms~^ for ry = 0.1 estimated by Favata et al. (2004), we 
estimate a recoil velocity of 100 it 20kms"^ for this mass ratio. For rj = 0.2 we estimate a recoil 
between 200 and 300kms~^, with a "best guess" of 250 km s^^ (the maximum velocity shown in 
Fig. 1 is 243kms~^). We regard our computation of the recoil in the adiabatic inspiral phase (up 
to the ISCO) as rather solid thanks to the accurate 2PN formula we use, and the fact that the 
1.5PN and 2PN results are so close to each other. However, obviously, using PN methods to study 
binary inspiral inside the ISCO is not without risks, and so it would be very desirable to see a 
check of our estimates using cither black hole perturbation theory (along the lines of Oohara & 
Nakamura (1983), Nakamura k, Haugan (1983) or Fitchett &: Detweiler (1984)) or full numerical 
relativity. It is relevant to point out that our estimates agree well with those obtained using 
numerical relativity in the "Lazarus approach" , or close- limit approximation, which treats the final 
merger of comparable-mass black holes using a hybrid method combining numerical relativity with 
perturbation theory (Campanelli 2005). In the small mass-ratio limit, they also agree well with a 
calculation of the recoil from the head-on plunge from infinity using perturbation theory (Nakamura 
et al. 1987). Therefore, wc hope that our estimates will enable a more focussed discussion of the 
astrophysical consequences of gravitational radiation recoil. 

The remainder of this paper provides details. In Section 2, we derive the 2PN accurate linear 
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momentum flux using a multipole decomposition, together with 2PN expressions for the multipole 

moments in terms of source variables. In Section 3 we speciaHze to binary systems, and to circular 
orbits. In Section 4, we use these results to estimate the recoil velocity and discuss various checks 
of our estimates. Section 5 makes concluding remarks. 



2. General formulae for linear momentum flux 

The flux of linear momentum P, carried away from general isolated sources, is first expressed 
in terms of symmetric and trace-free (STF) radiative multipole moments, which constitute very 
convenient sets of observables parametrizing the asymptotic wave form at the leading order |X|~^ 
in the distance to the source, in an appropriate radiative coordinate system = (T, X) (Thorne 
1980). Denoting by Ui^...ii,{T) and T^^...j^(T) the mass-type and current-type radiative moments at 
radiative coordinate time T (where I is the multipolar order), the linear momentum flux reads 



rm = Vf 2(£ + 2)(^ + 3) (1) ,rj.. 
•^P^^^ I ^(^ + 1)1(2^ + 3)!! ^ih-ky-' i^'h-ky-' ) 



+ ^(^ + ^) iT)V^^^ (T) 
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where the superscript (n) refers to the time-derivatives, and Sijk is Levi-Civita's antisymmetric 
symbol, such that £123 = Taking into account all terms up to relative 2PN order (in the case 
of slowly moving, PN sources), we obtain 

-ri _ 2 (1) (1) 16 (1) (1) 

1134 ^^''^ j''^ 126 63 ^■'^ 

59400 jKlm 14175 jlmn Klmn 945 y'''' J''-' 

The first two terms represent the leading order in the linear momentum flux, which corresponds 
to radiation reaction effects in the source's equations of motion occuring at the 3.5PN order with 
respect to the Newtonian force law. Indeed, recall that although the dominant radiation reaction 
force is at 2.5PN order, the total integrated radiation reaction force on the system (which gives 
the linear momentum loss or recoil) starts only at the next 3.5PN order (Peres 1962; Bonnor 
&; Rotenberg 1961; Papapetrou 1962). Radiation reaction terms at the 3.5PN level for compact 
binaries in general orbits have been computed by Iyer &; Will (1995), Jaranowski Sz Schafer (1997), 
Pati & Will (2002), Konigsdorffer et al. (2003) and Nissanke & Blanchet (2005). In Eq. (3) all the 
terms up to 2PN order relative to the leading linear momentum flux are included. This precision 
corresponds formally to radiation reaction effects up to 5.5PN order. 
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The radiative multipole moments, seen at (Minkowskian) future null infinity, Uij^...^ and 
Vi^...ii,, are now related to the source multipole moments, say and Jij^...i^, following the post- 

Minkowskian and matching approach of Blanchet & Damour (1986) and Blanchet (1995, 1998). 
The radiative moments differ from the source moments by non-linear multipole interactions. At 
the relative 2PN order considered in the present paper, the difference is only due to interactions of 
the mass monopole M of the source with higher moments, so-called gravitational-wave tail effects. 
For the source moments and Ji,^...j^.. wc use the expressions obtained in Blanchet (1995, 1998), 

valid for a general extended isolated PN source. These moments are the analogues of the multipole 
moments originally introduced by Epstein & Wagoner (1975) and generalized by Thorne (1980), 
and which constitute the building blocks of the direct integration of the retarded Einstein equations 
(DIRE) formalism (Will &; Wiseman 1996; Pati & Will 2000). The radiative moments appearing 
in Eq. (3) are given in terms of the source moments by (see Eqs. (4.35) in Blanchet (1995)) 
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where M = I denotes the constant mass monopole or total ADM mass of the source. The relative 
order of the tail integrals in Eqs. (4) is 1.5PN. The constant h entering the logarithmic kernel of 
the tail integrals represents an arbitrary scale which is defined by 



r = tH-pH-2Mln(^) , 



(5) 



where tn and pn correspond to a harmonic coordinate chart covering the local isolated source (pn is 
the distance of the source in harmonic coordinates). We insert Eqs. (4) into the linear momentum 
flux (3) and naturally decompose it into 
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where the "instantaneous" piece, which depends on the state of the source only at time T, is given 

by 
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and the "tail" piece, formally depending on the integrated past of the source, reads 
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The four terms in Eq. (8) correspond to the tail parts of the moments parametrizing the "Newto- 
nian" approximation to the flux given by the first line of (3). All of them will contribute at 1.5PN 
order. 



3. Application to compact binary systems 

We specialize the expressions given in Section 2, which are valid for general PN sources, to the 
case of compact binary systems modelled by two point masses mi and m2- For this application, 
all the required source multipole moments up to 2PN order admit known explicit expressions, 
computed in Blanchct ct al. (1995, 2002) and Arun et al. (2004) for circular binary orbits. Here we 
quote only the results. Mass parameters arc m = mi +m2, 5m = mi — m2 and the symmetric mass 
ratio rf = m\m2/m^. We define x = xi — X2 and r = |x| to be the relative vector and separation 
between the particles in harmonic coordinates, respectively, and v = dx./dt to be their relative 
velocity {t = tn is the harmonic coordinate time). We have, for mass- type moments. 
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and, for current-type moments. 
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We indicate the symmetric-trace-free projection using carets surrounding indices. Thus, the STF 
product of £ spatial vectors, say x'^'^""^^ = • • -x**, is denoted x^^^'"^''^ = STF [a:;*i"'*^]. Similarly, 
we pose = STF [x^i-ikyik+i-iey 

The total mass M in front of the tail integrals in (4) simply reduces, at the approximation 
considered in this paper, to the sum of the masses, i.e. M = m = nii + m2- Thus, to compute the 
tail contributions (8), we simply need the Newtonian approximation for all the moments. 

As seen in Eqs. (7)- (8) we need to perform repeated time-differentiations of the moments. 
These are consistently computed using for the replacement of accelerations the binary's 2PN equa- 
tions of motion in harmonic coordinates (for circular 2PN orbits) 

- = -a.V, (11) 

where u denotes the angular frequency of the circular motion, which is related to the orbital 
separation r by the generalized Kepler law 

-^ = 5{i + ?(-3 + .) + (7)X« + f'' + ''0}- <i^> 

The inverse of this law yields [using x = (muj)'^^^] 

= =x{l..(.-|)..= (l-f,)}. (13) 

The tail integrals of Eq. (8) are computed in the adiabatic approximation by substituting 
into the integrands the components of the moments calculated for exactly circular orbits, with the 
current value of the orbital frequency u (at time T), but with different phases corresponding to 
whether the moment is evaluated at the current time T or at the retarded time t < T. For exactly 
circular orbits the phase difference is simply Ac/) = u){T — t). All the contractions of indices are 
performed, and the result is obtained in the form of a sum of terms which can all be analytically 
computed by means of the mathematical formula 

where u is the orbital frequency, n the number of the considered harmonics of the signal (n = 1, 
2 or 3 at the present 2PN order) and C = 0.577- • • is Euler's constant. As shown in Blanchet & 
Schafer (1993) (see also Blanchet et al. (1995); Arun et al. (2004)), this procedure to compute the 
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tails is correct in the adiabatic limit, i.e. modulo the neglect of 2.5PN radiation reaction terms 
0[{v/c)^] which do not contribute at the present order. 

As it will turn out, the effect of tails in the linear momentum flux comes only from the first term 
in the right side of Eq. (14), proportional to tt. All the contributions due to the second term in (14), 
which involves the logarithm of frequency, can be reabsorbed into a convenient definition of the 
phase variable, and then shown to correspond to a very small phase modulation which is negligible 
at the present PN order. This possibility of introducing a new phase variable containing all the 
logarithms of frequency was usefully applied in previous computations of the binary's polarization 
waveforms (Blanchet et al. 1996; Arun et al. 2004). We introduce the phase variable ijj differing 
from the actual orbital phase angle 0, whose time derivative equals the orbital frequency = u), 

by 

= (f)-2mujln(^) , (15) 



where u) denotes a certain constant frequency scale that is related to the constant b which was 
introduced into the tail integrals (4), and parametrizes the coordinate transformation (5) between 
harmonic and radiative coordinates. The constants Q and b are in fact devoid of any physical 
meaning and can be chosen at will (Blanchet et al. 1996; Arun et al. 2004). To check this let us 
use the time dependence of the orbital phase </> due to radiation-reaction inspiral in the adiabatic 
limit, given at the lowest quadrupolar order by (see e.g. Blanchet et al. (1996)) 

where Tc and <pc denote the instant of coalescence and the value of the phase at that instant. Then 
it is easy to verify that an arbitrary rescaling of the constant a) by cD — *■ A tD simply corresponds to 
a constant shift in the value of the instant of coalescence, namely Tc —>■ Tc + 2m In A. Thus, any 
choice for Q is in fact irrelevant since it is equivalent to a choice of the origin of time in the wave 
zone. The relation between tD and b is given here for completeness, 
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The irrelevance of Q and b is also clear from Eq. (5) where one sees that they correspond to an 
adjustement of the time origin of radiative coordinates with respect to that of the source-rooted 
harmonic coordinates. 

Let us next point out that the phase modulation of the log-term in Eq. (15) represents in 
fact a very small effect, which is formally of order 4PN relative to the dominant radiation-reaction 
expression of the phase as a function of time, given by (16). This is clear from the fact that Eq. 
(16) is of the order of the inverse of radiation-reaction effects, which can be said to correspond 
to — 2.5PN order, and that, in comparison, the tail term is of order -I-1.5PN, which means 4PN 
relative order. In the present paper we shall neglect such 4PN effects and will therefore identify 
the phase tp with the actual orbital phase of the binary. 
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We introduce two unit vectors n* and A', respectively along the binary's separation, i.e. in the 
direction of the phase angle ip, and along the relative velocity, in the direction oi ip + ^, namely 



^ COSljj \ 



n = 







and A' 




(18) 



Finally, the reduction of the two terms (7) and (8) for compact binaries using the source moments 
(9)-(10) is straightforward, and yields the complete expression of the 2PN linear momentum flux, 
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The first term is the "Newtonian" one which, as we noted above, really corresponds to a 3.5PN 
radiation reaction effect. It is followed by the IPN relative correction, then the 1.5PN correction, 
proportional to tt and which is exclusively due to tails, and finally the 2PN correction term. We 
find that the IPN term is in agreement with the previous result by Wiseman (1992). The tail 
term at order 1.5PN and the 2PN term arc new with the present paper. Alternatively we can also 
express the flux in terms of the orbital frequency uj, with the help of the PN parameter defined by 
X = {mu)'^/^. Using Eq. (13) we obtain 
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The latter form is interesting because it remains invariant under a large class of gauge transforma- 
tions. 

Next, in order to obtain the local loss of linear momentum by the source, we apply the mo- 
mentum balance equation 

^ = -^p(r), (21) 

which yields Eq. (1). Upon integration, this yields the net change of linear momentum, say 
AP* = — dtTp{t). In the adiabatic limit, i.e. at any instant before the passage at the ISCO, 
the closed form of AP* can be simply obtained (for circular orbits) from the fact that dn^/dt = a; A* 
and the constancy of the orbital frequency u. This is of course correct modulo fractional error 
terms 0[{v/c)^] which are negligible here. So, integrating the balance equation (21) in the adiabatic 
approximation simply amounts to replacing the unit vector A* by W and dividing by the orbital 
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frequency u. In this way we obtain the recoil velocity = IS.P'^ jm as' 
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or, alternatively, in terms of the x-parameter. 
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Equations (1) and (23) will be the basis for our numerical estimates of the recoil velocity, to be 
carried out in the next Section. 



4. Estimating the recoil velocity 

4.1. Basic assumptions and analytic formulae 

We now wish to use Eqs. (1) and (23) to estimate the recoil velocity that results from the 
inspiral and merger of two black holes. It is clear that the PN approximation becomes less reliable 
inside the innermost stable circular orbit (ISCO). Nevertheless, we have an expression that is 
accurate to 2PN order beyond the leading effect, which will therefore be very accurate over all the 
inspiral phase all the way down to the ISCO, so we have some hope that, if the higher-order terms 
can be seen to be small corrections throughout the process, we can make a robust estimate of the 
overall kick. 

In Eq. (23) we have re-expressed the recoil velocity in terms of the orbital angular velocity 
CO, Eq. (12), consistently to 2PN order. One advantage of this change of variables is that the 
momentum loss is now expressed in terms of a somewhat less coordinate dependent quantity, namely 
the orbital angular velocity as seen from infinity. A second advantage is that the convergence of 

the PN scries is significantly improved. In terms of the variable m/r, the coefficients of the IPN 
and 2PN terms arc of order —10 and 33 - 41, respectively, depending on the value of ry, whereas in 
terms of x, they are of order —5 and —3 — 1-1.4, respectively. 

We assume that the system undergoes an adiabatic inspiral along a sequence of circular orbits 
up to the ISCO. For the present discussion the ISCO is taken to be the one for point-mass motion 



''The recoil could also be defined from the special-relativistic relation = AP^ / \''m? + AP-^, but since AP* is of 
order 3.5PN the latter "relativistic" definition yields the same 2PN results, and in fact differs from our own definition 
by extremely small corrections, at the 7PN order. 
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around a Schwarzschild black hole of mass m = mi + m2, namely mwisco = 6 or xisco = V^- 
The recoil velocity at the IS CO is thus given by 



ISCO 



464 Sm 2 4 
105—^ ^ISCO 



1 + 



1139 



+ 



71345 



+ 



452 

36761 147101 



309 



3/2 



22968 2088 



-1] + 



68904 



2 \ 2 
^ ) ^ISCO 



"isco 



(24) 



In order to determine the kick velocity accumulated during the plunge, we make a number of 
simplifying assumptions. We first assume that the plunge can be viewed as that of a "test" particle 
of mass n moving in the fixed Schwarzschild geometry of a body of mass m, following the "effective 
one-body" approach of Buonanno & Damour (1999) and Damour (2001). We also assume that the 
effect on the plunge orbit of the radiation of energy and angular momentum may be ignored; over 
the small number of orbits that make up the plunge, this seems like a reasonable approximation 
(Favata et al. (2004) make the same assumption). 



We therefore adopt the geodesic equations for the Schwarzschild geometry. 



dt 
J? 



E 



1 — 2m/rs ' 
L 



drs 
dr 



E' 



2m 
rs 



1 + 



L2 



(25a) 
(25b) 
(25c) 



where r is proper time along the geodesic, E is the energy per unit mass (fi in this case), and 
L = mL is the angular momentum per unit mass. Then, from Eqs. (25b) and (25c), we obtain the 
phase angle of the orbit as a function of y = m/rs by 



1/2 



dy, 



,y, v£;2_(i_2y)(i + LV)- 
where we choose -0 = at the beginning of the plunge orbit defined by y = yo- 
The kick velocity accumulated during the plunge is then given by^ 



(26) 



plunge 



— / —dt. 

m J to at 



(27) 



However, the coordinate time t is singular at the event horizon, so we must find a non-singular 
variable to carry out the integration. We choose the "proper" angular frequency. Hi = dip /dr. In 



''The radiative time T in the linear momentum loss law (21) can be viewed as a dummy variable, and we henceforth 
replace it by the Schwarzschild coordinate time t. 
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addition to being monotonically increasing, this variable has the following useful properties along 
the plunge geodesic: 

mu = Ly^ , (28a) 
muj = mu)^—^ = ^y^{l-2y), (28b) 

^ = l^y\E^-{l-2y){l+-LW)X" . (28c) 
at m I J 



Then 



^ plunge 



1 f dP' du! 



m J dt dijj/dt 

/"l/Horizon / ]^ 



dP" \ dy 

* ^ (29) 



£;2-(l-2y)(l + LV: 



1/2 



where yo is defined by the matching to a circular orbit at the ISCO that we shall discuss below. 

Notice that, because dP^/dt oc x^^/^ oc (mo;)^^/^, the quantity in parentheses in Eq. (29) is 
well behaved at the horizon; in fact it vanishes at the horizon because a; = there [cf. Eq. (28b)]. 
Thus, we find that the integrand of Eq. (29) behaves like (ma;)^/^ oc (1 — 2y)^/^ at the horizon, 
and the integral is perfectly convergent. Furthermore, since the expansion of dP^/dt is in powers of 
muj, the convergence of the PN series is actually improved as the particle approaches the horizon. 
To carry out the integral, then, we substitute for x = (mco)^^^ in dP^ /dt using Eq. (28b), and 
integrate over y. 

We regard this approach as robust, because it uses invariant quantities such as angular fre- 
quencies, and uses the nature of the flux formula itself to obtain an integral that is automatically 
convergent. Favata et al. (2004) tried to control the singular behavior of the t integration with an 
ad hoc regularization scheme. 

We then combine Eqs. (24) and (29) vectorially to obtain the net kick velocity, 

Al^^ = Vi^co + AFpV„g,, (30) 

in which V^^qq is given by Eq. (24) above with nJgQQ = (1,0,0). 

There are many ways to match a circular orbit at the ISCO to a suitable plunge orbit; we use 
two different methods. In one, we give the particle an energy E such that, at the ISCO, and for an 
ISCO angular momentum Zisco = \/l2m, the particle has a radial velocity given by the standard 

quadrupole energy-loss formula for a circular orbit, namely dr^/dt = — (64/5)?7(m/rH)^, where rn 
is the orbital separation in harmonic coordinates. At the ISCO for a test body, rn = 5m, so we 
have {dr}i/dt)isco = -{8/25)^7]. This means also {drs/dt)isco = "(8/25)^77 in the Schwarzschild 
coordinate rs = rn + m (recall that ts = tu = t). It is straightforward to show that the required 
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energy for such an orbit is given by 



9 fdrsV 



4:\dt J IgQQ 



(31) 



We therefore integrate Eq. (29) with that energy, together with Lisco = vl2 and the initial 
condition yo = yisco = 1/6 (from Eq. (28b) we note that, with this choice of initial condition, 
muQ ^ 6""^/^). We choose also to terminate the integration when rs = 2(m + /x) hence ^Horizon ~ 
2(1 + r/). 

With this initial condition, the number of orbits ranges from 1.2 for 77 = 1/4 to 1.8 for 77 = 1/10 
to 4.3 for ?7 = 1/100. It is also useful to note that the radial velocity remains small compared to 
the tangential velocity throughout most of the plunge; the ratio {dro, / dr) / {ro, dijj/d.T) = r^^drs/dtl^ 
reaches 0.14 at rs = 4m, 0.3 at rs = 3m, and 0.5 at rs = 2m, roughly independently of the 
value of 77. This justifies our use of circular orbit formulae for the momentum flux as a reasonable 
approximation. 

In a second method, we evolve an orbit at the ISCO piecewise to a new orbit inside the ISCO, 
as follows: using the energy and angular momentum balance equations for circular orbits in the 
adiabatic limit at the ISCO, we have 

dE 32 r] . , , 

-dF = -Ym^^^cO' ^''"^ 

^ - ^--^-^nx'/' (32b) 

We approximate these relations by "discretizing" the variations of the energy and angular momen- 
tum in the left sides around the ISCO values Ejsco = v^8/9 and I/isco = V^m. Hence, we write 
dE/dt = {E — Eisco)/{(xP) and dL/dt = {L — Ljsco) / {ctP) , where aP denotes a fraction of the 
orbital period P of the circular motion at the ISCO. Using then wisco = 27r/P = 'm^^x^^QQ this 
gives the following values for the plunge orbit 

E = Eisco - ^ a V a^isco ' (3^*) 
L = -^isco - ^ a rj xfsco ■ (33b) 

Then, in this second model we integrate Eq. (29) with the latter values, and using the initial 
inverse radius yo = ('Ti/rs) initial of this new orbit which is given by the solution of the equation 

mwisco = 6-^/2 = i y^il - 2yo) . (34) 
E 

For the final value W6 simply take tlic liorizon at rg — 2?77. (licnce ^/Horizon 

= 1/2), in the spirit 

of the effective one-body approach (Buonanno & Damour 1999; Damour 2001) where the binary's 
total mass m is identified with the black-hole mass and where /j, is the test particle's mass. For the 
fraction a of the period, we choose values between 1 and 0.01, and check the dependence of the 
result on this choice (see below). 
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Table 1: Recoil velocity (kms ) at the ISCO defined by xisco = 1/6. 



77 = pL/m 


0.05 


0.1 


0.15 


0.2 


0.24 


Newtonian 


2.29 


7.92 


14.56 


18.30 


11.78 


N + IPN 


0.27 


0.77 


1.16 


1.12 


0.55 


N + IPN + 1.5PN (tail) 


2.87 


9.80 


17.74 


21.96 


13.97 


N + IPN + 1.5PN + 2PN 


2.73 


9.51 


17.57 


22.22 


14.38 



4.2. Numerical results and checks 

First, we display the recoil velocities at the ISCO given by Eq. (24) for each PN order and 
various values of 77 in Table 1. The 2PN values of the velocity at the ISCO are also plotted as a 
function of 77 in Figure 1 (dot-dash curve). On should note, from Table 1, the somewhat strange 
behavior of the IPN order, which nearly cancels out the Newtonian approximation (as already 
pointed out by Wiseman (1992)). The maximum velocity accumulated in the inspiral phase is 
around 22kms^^. 

Next, we evaluate the kick velocity from the plunge phase, and carry out a number of tests of 
the result. In our first model, where the plunge energy is given by (31), we choose rg = 6m as the 
ISCO, and rg = 2(m + /i) = 27ti(1 + rf) as the final merger point. The latter value corresponds 
to the sum of the event horizons of black holes of mass m and //, and is an effort to estimate the 
end of the merger when a common event horizon envelops the two black holes, and any momentum 
radiation shuts off. 

The resulting total kick velocity as a function of rj is plotted as the solid (red) curve in 
Figure 1. We also consider the kick velocity generated when we take only the leading "Newtonian" 
contribution (dashed [black] curve), and when we include the IPN terms (short dashed [green] 
curve) and the IPN + 1.5PN terms (dotted [blue] curve). Notice that, because the IPN term has a 
negative coefficient, the net kick velocity at IPN order is smaller than at Newtonian order. On the 
other hand, because the 2PN coefficient is so small, the 1.5PN correct value and the 2PN correct 
value arc very close to each other. 

In order to test the sensitivity of the result to the PN expansion, we have considered terms 
of 2.5PN, 3PN and 3.5PN order, by adding to the expression (1) terms of the form a2.5PNa^^''^ + 
asPN^c^ + os.sFN^c^'^^) and varying each coefficient between +10 and —10. For example, varying 
a2.5PN and ospn, leads to a max;imum variation in the velocity of ±30% [i.e. between the values 
(—10,-10) and (10,10)] for a range of 77. Assuming that the probability of occurrence of a specific 
value of each coefficient is uniform within the interval [—10,10], we estimate an rms error in the 
kick velocity, shown as "error bars" in Figure 1. Varying aa.sPN between —10 and 10 has only a 
10% effect on the final velocity. These considerations lead us to crudely estimate that our results 
are probably good to ±20%. 
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In the limit of small 77, our numerical results give an estimate for the kick velocity: 

— « 0.043 v^l -4r? when r/ ^ , (35) 

with the coefficient probably good to about 20 %. 

We also test the sensitivity of the results to the end point: carrying out the integration all 

the way to rs = 2m, as in our second model, Eqs. (32)-(34), has only a one percent effect on the 
velocity for r\ = 0.2, and has essentially negligible effect for smaller values of rj. Wc also vary the 
value of the radius where we match the adiabatic part of the velocity with the beginning of the 
plunge integration. For matching radii between 5.3 m and 6 m, the final kick velocity varies by at 
most seven percent for 77 = 0.2 and five percent for 77 = 0.1. 

In establishing the initial energy for the plunge orbit, we used the quadrupole approximation 
for drn/dt in harmonic coordinates. We have repeated the computation using a 2PN expression 
for drji/dt expressed in terms of mw; the effect of the change is negligible. 

Our second method for matching to the plunge orbit, Eqs. (32)-(34), gives virtually identical 
results. For the 2PN correct values, and for values of the parameter a below 0.1, this method gives 
velocities that are in close agreement with those shown in Figure 1. For instance, with a = 0.1 and 
r] = 0.2, the kick velocity is equal to 245kms~^, compared to 243kms~^ with the first method. 
Small values of a correspond to a smoother match between the circular orbit at the ISCO and the 
plunge orbit. For a = 1, implying a cruder match, the kick velocities are lower than those shown 
in Figure 1: 4 % lower for 77 = 0.1, 10 % lower for rj = 0.2, and 14 % lower for 77 = 0.24. These 
differences are still within our overall error estimate of about 20 % indicated in Figure 1. 



5. Concluding remarks 

Our results are consistent with, but substantially sharper than the estimates for kick velocity 
for non-spinning binary black holes given by Favata et al. (2004). They are also consistent with 
estimates given by Campanelli (2005) obtained from the Lazarus program for studying binary black 
hole inspiral using a mixture of perturbation theory and numerical relativity. A recent improved 
analysis (Campanelh 2005) gives 240 ± 140 km s"^ at 77 = 0.22 and 190 ± 100 km s"^ at rj = 0.23; as 
compared with our estimates of 211ib40 and 183 it 37, respectively. In the limit of small mass ratio, 
Eq. (35) agrees very well with the result V/c = 0.04577^ obtained by Nakamura et al. (1987) using 
black hole perturbation theory for a head-on collision from infinity. Since, as we have seen, the 
contribution of the inspiral phase is small and the recoil is dominated by the final plunge, one might 
expect a calculation of the recoil from a head-on plunge to be roughly consistent with that from a 
plunge following an inspiral, despite the different initial conditions; accordingly the agreement we 
find with Nakamura et al. (1987) for the recoil values is satisfying^. 



®We thank Marc Favata for bringing the Nakamura et al. (1987) work to our attention 
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Finally, we remark on the curious fact that our 2PN result shown in Figure 1 can be fit to 
better than one percent accuracy over the entire range of 77 by the simple formula 



While we ascribe no special physical significance to this formula in view of the uncertainties in our 
PN expansion, it illustrates that, beyond the overall rf y/\ — Ar] dependence, the post-Newtonian 
corrections and the plunge orbit generate relatively weak dependence on the mass ratio. Such an 
analytic formula may be useful in astrophysical modeling involving populations of binary black hole 
systems. 

Inclusion of the effects of spin will alter the result in several ways. First, it will allow a net 
kick velocity even for equal mass black holes. Second, it will significantly change the plunge orbits, 
depending on whether the smaller particle orbits the rotating black hole in a prograde or retrograde 
sense. In futTirc work, wc plan to treat this problem using our 2PN formulae for linear momentum 
flux, augmented by the 1.5PN spin orbit flux terms of Kidder (1995), combined with a similar 
treatment of plunge orbits in the equatorial plane of the Kerr geometry. 

We are grateful to Jeremiah P. Ostriker for motivating us to look at this question and for useful 
discussions, to Masaru Shibata for useful input related to black hole perturbation calculations of 
momentum recoil, and to Scott Hughes and Marc Favata for useful comments on the manuscript. 
This work is supported in part by the US National Science Foundation under grant No. PHY 
03-53180. CW is grateful for the support of the Centre National dc la Recherche Scientifique, and 
the hospitality of the Institut d'Astrophysique de Paris, where this work was completed. MSSQ 
thanks B. R. Iyer and K. G. Arun for useful discussions. The calculations were performed with the 
help of the softwares Mathematica and Maple. 
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